Mechanical properties of polycrystalline graphene based on a realistic atomistic model 
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Graphene can at present be grown at large quantities only by the chemical vapor deposition 
method, which produces polycrystalline samples. Here, we describe a method for constructing 
realistic polycrystalline graphene samples for atomistic simulations, and apply it for studying their 
mechanical properties. We show that cracks initiate at points where grain boundaries meet and then 
propagate through grains predominantly in zigzag or armchair directions, in agreement with recent 
experimental work. Contrary to earlier theoretical predictions, we observe normally distributed 
intrinsic strength (^ 50% of that of the mono-crystalline graphene) and failure strain which do not 
depend on the misorientation angles between the grains. Extrapolating for grain sizes above 15 nm 
results in a failure strain of ^ 0.09 and a Young's modulus of ^ 600 GPa. The decreased strength 
can be adequately explained with a conventional continuum model when the grain boundary meeting 
points are identified as Griffith cracks. 
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Grain boundaries define the electronic and mechani- 
cal properties of polycrystalline materials. Since chemi- 
cal vapor deposition (CVD) is currently the only way for 
producing industry-scale graphene membranes, and leads 
to polycrystalline samples, study of grain boundaries in 
graphene has become of fundamental importance during 
the recent years. In a two-dimensional material, such 
as graphene, the boundaries also have a critical contri- 
bution to the chemical reactivity. Because of this, al- 
though atomic scale imaging can in principle reveal their 
exact structure, the boundaries tend to be covered by 
adsorbates with only short segments available for direct 
imaging. Nevertheless, experiments [1-6] have revealed 
meandering serpent-like boundaries which are typically 
formed from pentagon-heptagon-pairs in the parts not 
covered by the adsorbates. 

Mechanical properties of graphene sheets have been a 
topic of intense research already for two decades in the 
context of carbon nanotubes (see Ref. for a topical 
review). More recently, in 2007 [8'], Liu and co-workers 
utilized ah initio calculations to study the elastic moduli 
and fracture characteristic of graphene. Young's mod- 
ulus was found to be 1.05 TPa, and failure strain, de- 
pending on the pulling direction, 0.194-0.266. Intrinsic 
strength was estimated to be 110-121 GPa, similarly de- 
pending on the pulling direction. The role of pre-existing 
defects on these properties was also studied [9[. It was 
noticed that their effect does not depend on the exact 
atomic structure of the defects but rather on their size. 
The authors also showed that the intrinsic strength of 
graphene with crack-like defects can be described with 
a continuum model using the Griffith formula for defect 
sizes down to 10 A. Soon after this, Frank et al. used 
a tip of an atomic force microscope to obtain a Young's 
modulus of 0.5 TPa for suspended stacks of graphene 
sheets [10]. A year later, Lee and co-workers reported on 
several mechanical properties of graphene using a similar 
technique [ll| , establishing graphene as the strongest ma- 
terial ever measured. They reported an intrinsic strength 



of 42 N/m (corresponding to 130 GPa assuming graphene 
thickness to be the inter-layer distance in graphite, i.e., 
0.335 nm) occurring at 0.25 strain. Young's modulus was 
estimated to be 1 TPa, in a good agreement with the- 
ory [8]. In 2009, Xiao and others reported on their the- 
oretical work p^], in which they obtained failure strain 
of ca. 0.10 for graphene sheets with Stone- Wales defects 
(one rotated bond) with an intrinsic strength very close 
to that of the pristine structure (the difference was larger 
for small-diameter nanotubes). This result would be con- 
sistent with the continuum model [9] assuming the defect 
corresponds to a crack with a size below 5 A, which seems 
reasonable for this defect. 

These early works concentrated on either mechan- 
ically exfoliated pristine graphene or mono-crystalline 
graphene with point defects. The first experimental stud- 
ies on mechanical properties of polycrystalline graphene 
samples were carried out only last year (2011) [2, 13, ij. 
The experimental results can be summarized as fol- 
lows: The intrinsic strength for polycrystalline samples 
is somewhat above one third of that for mono-crystalline 
graphene (ca. 35 GPa) [2|, 113| and cracks propagate 
through the bulk of the grains (l4J mostly along zigzag 
and armchair directions, not along grain boundaries as 
could be intuitively expected. 

In the meanwhile, also theoretical work on energet- 
ics and other non-elastic properties [5|, Il5l - [l7j of grain 
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boundaries as well as on their mechanical response 13, 
18] has been carried out. Total energy calculations [16|, 
17] have established that the idealized low-energy con- 
figuration of grain boundaries is a linear tilt boundary 
consisting of a repeating set of pent agon- heptagon-pairs 
which act as dislocation cores in a lattice otherwise con- 
structed of hexagonal carbon rings. Intrinsic strength of 
graphene sheets with infinitely long such boundaries has 
been estimated to be 50-100 GPa with failure strains in 
the range of 0.07-0.15 [18], depending on the misorien- 
tation angle between the adjacent grains. Since a higher 
misorientation angle yields a higher density of disloca- 
tion cores at the linear tilt boundary, but also to higher 
intrinsic strength of the model structures, the authors 
noted that their results disagree with the fracture me- 



chanics model, assuming the heptagons of the dislocation 
cores correspond to Griffith cracks, which would predict 
graphene sheets to become weaker with an increasing de- 
fect density. If the grain boundaries themselves are in- 
deed the weakest point in the lattice, this can be argued 
to be a reasonable comparison since also Stone- Wales de- 
fect consists of pentagons and heptagons and it has been 
shown to weaken graphene [12]. However, despite simi- 
larities between the theoretical models [16|, |T3 and short 
segments of the actually observed non-decorated bound- 
aries [l|-[5| , it remains unclear whether such infinitely long 
linear arrangements of dislocation cores can serve as a re- 
alistic model for studying mechanical properties of poly- 
crystalline graphene. 

Theoretical calculations presented along the experi- 
mental work in Ref. jl3| assumed that voids would exist 
in polycrystalline graphene samples and that they could 
explain the apparent discrepancy between the experimen- 
tal results of ~ 35 GPa and the theoretical estimates of 
50-100 GPa. However, it is questionable how well this 
model corresponds to actual polycrystalline graphene 
samples. Moreover, the inherent difficulties in assessing 
the mechanical properties of a membrane suspended on 
a hole by applying force with a tip of a microscope ne- 
cessitate theoretical confirmation with a realistic model 
system. 

Here, we describe an automated method for creating 
polycrystalline graphene structures with realistic misori- 
entation angle and carbon ring size statistics as well as 
serpent-like boundaries similar to those observed exper- 
imentally. Using atomistic simulations, we then sub- 
ject our samples to a study of their mechanical prop- 
erties. We show that close to the failure strain, cracks 
appear typically at the points where grain boundaries 
meet, and in agreement with the recent experimental 
studies, then propagate through grains predominantly in 
zigzag or armchair directions. Contrary to earlier theo- 
retical predictions fl8] , neither intrinsic strength nor fail- 
ure strain of our samples depend on the misorientation 
angle between the grains, but are normally distributed 
similar to recent experimental studies, [2|, [l3| where in- 
trinsic strength of ~ 35 GPa was reported. We obtain 
a slightly higher value (46 GPa) which is still in a rea- 
sonably good agreement with the experimental one. At 
the large grain size limit (^ 15 nm) the failure strain 
is about 0.09 and Young's modulus is close to 600 GPa. 
The formation of cracks at the meeting points of grain 
boundaries, completely neglected in the previous theo- 
retical studies, resolves the discrepancy between the ex- 
periments and the theoretical results and shows that the 
Griffith model can after all be used to describe the me- 
chanical properties of polycrystalline graphene samples 
when a realistic atomistic model is used. 

Without pre-patterned seeds for growth, CVD growth 
of graphene is initiated at several nucleation sites simul- 
taneously. On a substrate such as Cu, which doesn't offer 
epitaxiality, the lattice orientations of the growing grains 
are random [^. To mimic such growth, we first wrote 



a computer code which creates a pre-selected number of 
randomly placed nucleation sites on a plane with pre- 
defined dimensions. For each such nucleation site {i) a 
random orientation Oi is selected for the lattice. In order 
to obtain approximately uniform size distribution for the 
grains, the sites are required to be at least 5.0 A apart 
from each other (5 A was selected arbitrarily). Next, we 
carry out an iterative process in which any of the missing 
neighbors of the already inserted atoms can appear with 
the same probability. When two grains approach, we use 
the following condition for deciding whether a lattice site 
is available for another atom: if d < 1.0 A or A^ > 3 {d 
is the distance between the lattice site and the closest 
existing atom, N is the number of atoms created closer 
than a — 1.0 A to the present site, where a is the length of 
the graphene lattice vector) the site is not free and will 
thus be disregarded for further growth. Upon testing, 
this condition was found to minimize the dangling bond 
density at the boundaries. 
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Figure 1. (Color online) Model structures for polycrystalline 
graphene. (a) Top and (b) side view of a periodic 20 nm 
X 20 nm graphene sheet with four grains, as marked by the 
numbered shaded areas. The lines indicate orientations of 
the graphene lattice within each grain. Note that all the pre- 
sented mechanical studies have been carried out for bicrys- 
talline samples, (c) Distribution of misorientation angles for 
the bicrystalline sample structures used in this study, (d) 
Relative probabilities for non-hexagonal carbon rings in the 
same structures. 

To equilibrate thus created polycrystalline samples, we 
first annealed them at 3000 K for 50 ps after which the 
system was quenched during a 10 ps run allowing the 
lattice to obtain its equilibrium size (pressure driven to 
zero). At this point the lattice appears somewhat crum- 
pled even after pressure relaxation since we did not re- 
strict relaxation in the out-of-plane direction. All sim- 
ulations were carried out with the classical molecular 
dynamics (MD) code PARCAS pQl - [22| with a reactive 
bond order potential developed by Brenner et at. [23|. 
Due to the large number of atoms in the structures 



(up to almost 10000) and large number of structures 
(385 in total), this is the only feasible method for car- 
rying out the simulations. A similar simulation setup 
has been used in earlier theoretical studies of mechani- 
cal properties of graphene [18], where a good agreement 
with ab initio methods has been noticed. Temperature 
and pressure control were handled using the Berendsen 
method [24]. The equilibration procedure leads to grain 
boundary structures similar to the 20 nm x 20 nm model 
presented in Fig. [TJi,b. 

After establishing a method for creating model struc- 
tures for polycrystalline graphene, we applied it for creat- 
ing 385 bicrystalline structures with grain sizes between 
~ 3-16 nm. Two randomly placed seeds were used for 
each structure to obtain exactly one misorientation angle 
(0) per structure. The resulting distribution of is pre- 
sented in Fig. [It, where i9 = 6>' = |6>i - 6>2| if (9' < 30° and 
(9 = 60° - 6>' otherwise (for graphene any Oi G [0°,60°]). 
As expected for two randomly selected orientations, the 
distribution is uniform with fluctuations resulting only 
from the finite sample size. In Fig. [l]i, we plot the 
relative occurrence of carbon rings other than hexagons 
within the created structures. The combined likelihood 
for tetragons and pentagons is similar to that of hep- 
tagons and octagons indicating mostly saturated bonds 
at the boundaries. The significantly lower probability 
for rings with more than eight atoms is a sign of an ex- 
isting but small local density deficit at the boundaries. 
Overall, the ring statistics seem reasonable. We noticed 
only very rarely if ever four-coordinated atoms which 
would indicate problems with the interaction model (such 
coordination is never observed in sp^-honded graphene 
even when it is heavily amorphized under an electron 
beam [25|, [26|). However, almost all structures contain a 
few under-coordinated atoms which could serve as reac- 
tive sites for covalently bonding adsorbates on the grain 
boundaries. Thus, the experimentally observed high cov- 
erage of grain boundaries gives further credibility for our 
model structures. 

Next, we subjected the created sample structures to 
extensive tensile testing. The simulations were carried 
out at 300 K as follows: we applied uni- axial strain in 
a step-by-step fashion always equilibrating the structure 
for 5 ps before increasing the strain (much slower pulling 
was also tested with no apparent change in the results). 
We employed periodic boundary conditions for all simu- 
lations in the in-plane {x and y) directions. For this part 
of the simulations we modified the cutoff of the interac- 
tion model to remove the unphysical softening at longer 
inter-atomic distances (above 1.92 A), which is crucial 
for many MD simulations, but in the present case only 
affects at high strains by leading to nonphysical features 
in the stress-strain curve. For inter-atomic distances be- 
low 1.92 A, the interaction model remained unchanged. 
Modifying the cutoff has also before been noticed to be 
required for properly addressing the mechanical proper- 
ties of graphene close to the fracture strain [18]. The con- 
tinuity of atomic trajectories and conservation of energy 
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Figure 2. (Color online) Stress-strain curves for polycrys- 
talline graphene samples as plotted for grain sizes of (a) 
d ^ 12 nm, (b) d ~ 6 nm and (c) d ~ 3 nm. The size of 
the grains is determined as an average diameter for a grain 
assumed to be circular as indicted in panel (d). 



were monitored during the simulations to avoid any prob- 
lems resulting from this modification. We also carried out 
test simulations to check that the differences caused to 
the stress-strain curves - either in pristine graphene or 
our bicrystalline samples - were limited to the unphysical 
features near fracture. 

The stress during deformation was calculated from the 
Virial expression as explained in Ref. [18] assuming a 
thickness of 0.335 nm for the graphene membrane. The 
resulting stress-strain curves for grains with three differ- 
ent average sizes are presented in Fig. [2l For pristine 
graphene, our method yields an intrinsic strength of 90- 
100 GPa at a failure strain of 0.15-0.20, depending on 
the pulling angle, in a good agreement with the ab ini- 
tio results [8]. As can be seen in Fig. [21 the presence 
of grain boundaries leads to approximately a 50% reduc- 
tion of the strength of the material independent of the 
grain size (intrinsic strength corresponds to the maxi- 
mum stress before the failure). For grain sizes above 
~ 12 nm, the average fracture strain is close to 0.10 
(Fig. [2^), whereas for the smaller ones it increases grad- 
ually (see Figs.[2]3,c) up to ca. 0.15 for the smallest rea- 
sonable grain sizes (~ 3 nm). This is because the grain 
boundaries are more fiexible than the bulk of the grains, 
and their role is pronounced at small grain sizes allow- 
ing higher overall strains. We point out that the Poisson 
effect has not been taken into account in the presented 
data. However, we checked whether it would affect the 
results by carrying out a subset of the simulations also 
without periodic boundaries in the ^/-direction. The ob- 
served deviations were within the uncertainties stemming 
from the finite sample size (that is, those seen in Fig. [2]). 

To better understand how the different measured prop- 
erties depend on the grain size {d) and misorientation be- 
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Figure 3. (Color online) Mechanical properties of the sample 
structures as a function of the misorientation between the 
two grains and the grain size d (or system size), (a) Failure 
strain, (b) intrinsic strength and (c) Young's modulus as a 
function of for d ^ 12 nm grains, (d) Failure strain as a 
function of d. (e) Distribution of the intrinsic strengths for 
all d. (f) Young's modulus as a function of d. 



tween the grains (^), we plot in Fig. [3] the failure strain, 
intrinsic strength and Young's modulus as functions of 
and d. For the data plotted as a function of (Fig. [3] 
a-c), we used only one grain size {d ^ 12 nm) to ease the 
interpretation of the results. What can be readily ob- 
served is that none of the calculated properties depend 
on 0. Instead, they are normally distributed over an av- 
erage value, in a stark contrast to the earlier theoretical 
prediction based on infinitely long linear grain boundary 
structures [18]. However, when the failure strain is plot- 
ted as a function of d (Fig. [3]i), a clear size-dependency 
emerges, as was qualitatively described above. For intrin- 
sic strength (Fig. [3^) we observe no d-dependency at all. 
Instead, the data is normally distributed for all d around 
a value of ~ 46 GPa. While Young's modulus is defined 
as the change in stress divided by the change in strain for 
the linear part of the stress-strain curve (Fig. [2j), its d- 



dependency is determined by that of the strain (Fig.[3f). 
For the strain, we can describe the different contributions 
of the bulk of the grain and the grain boundary with a 
constant describing the large-grain-size-limit eo and term 
inversely proportional to the grain size {oc d~^) to obtain 
e = eo + const. /d. Fitting this equation to the data for 
the failure strain yields a very good agreement as can be 
seen in Fig. [3li. Through the fit we obtain failure strain 
for large grains of ^ 0.09. A similar fit for the Young's 
modulus (Fig. [3]F) gives a value of ~ 600 GPa similarly 
for large grains. 

To further understand why the breaking stress is d- 
independent, we visually analysed the evolution of the 
atomic structure of our samples upon fracture. An ex- 
ample is presented in Fig. |4j What we noticed is that 
the crack formation often occurs at the points where the 
grain boundaries meet (marked with a square in Fig.|4^). 
After the crack is formed, however, it propagates typi- 
cally along the armchair or zigzag lattice directions in 
the bulk of the neighboring grains, similar to what has 
been recently suggested based on experimental observa- 
tions [J]. While the atomic structure of the boundaries, 
and that of their meeting points, are independent on the 
grain size, the fracture properties must also be grain size- 
independent, which is exactly what we observe in our re- 
sults. Moreover, the characteristic size of these meeting 
points in our structures (as can be seen in Fig. [TJi and 
Fig. [4^) is ~ 2 nm, which can be compared to the Griffith 
model data from Ref. [9]. The data would indicate that 
such a crack size would result in roughly 50% reduction 
in the strength of the material, in agreement with our 
data. 

As a conclusion, we established a method for creat- 
ing realistic polycrystalline graphene samples for atom- 
istic simulations. We applied this method for creating a 
representative set of samples for mechanical testing, and 
showed that the presence of grain boundaries reduces the 
strength of graphene by about 50% (down to ~ 46 GPa), 
in a reasonable agreement with experiments [2|, |l3| . How- 
ever, we observed no misorientation-dependency on any 
of the mechanical properties of the created samples which 
was recently suggested based on a theoretical study 
on graphene structures with infinitely long linear grain 
boundaries [l8|. Furthermore, we showed that crack for- 
mation occurs at points at which the boundaries meet 
and that the cracks propagate through the bulk of the 
neighboring grains typically along armchair and zigzag 
directions, similar to recent experimental findings [^]. 
The failure strain for polycrystalline graphene with grain 
sizes :^ 15 nm was found to be ~ 0.09 with a correspond- 
ing Young's modulus of ~ 600 GPa. Overall our results 
show that the mechanical properties of polycrystalline 
graphene can be reasonably well described using the con- 
tinuum model if the grain boundary meeting points are 
identified as the Griffith cracks in this material. 
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Figure 4. (Color online) An example case of crack formation and propagation, (a) The structure with no strain with different 
grains (1 and 2) marked with different colors. The lines indicate the orientations of the grains. The square shows the area 
where the grain boundaries meet and the crack will appear, (b-e) Snapshots of the area immediately around the crack during 
straining showing how the crack penetrates along the zigzag axis of the bulk of the grain, (f) The structure after the crack has 
penetrated through the grains. 
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